Mapping fire damage from the wildfires in Larache Province, Morocco, in 2022 using Sentinel-2 Normalised Burn Ratio (NBR) and differenced NBR1

Connecting to Google Drive from Colab.

# Loading the Drive helper and mount your Google Drive as a drive in the virtual machine
from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

Importing all required libraries

# Installing some libraries that are not on Colab by default
!pip install rasterio
!pip install geopandas
!pip install rasterstats
!pip install earthengine-api
!pip install requests
!pip install sentinelsat

# Importing libraries
import geopandas as gpd
import rasterio
from rasterio import plot
from rasterio.plot import show_hist
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal, ogr, osr
import json
import os
from os import listdir
from os.path import isfile, isdir, join
import math
from pprint import pprint
import shutil
import sys
import zipfile
import requests
import io
import webbrowser
import ee
import pandas as pd
import rasterio.mask
import xarray as xr
import matplotlib.colors as colors
import matplotlib
import glob
from rasterio.plot import show
from rasterio.mask import mask
from shapely.geometry import mapping

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# defining the root directory where our data are to be stored 
rootdir = '/content/drive/MyDrive/satelliteCW1' # this is where pygge.py is saved on my Google Drive

if rootdir not in sys.path:
    sys.path.append(rootdir)
# importing the pygge library of functions for this module
import pygge

%matplotlib inline

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting rasterio
  Downloading rasterio-1.3.6-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (20.0 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 20.0/20.0 MB 104.8 MB/s eta 0:00:00
Collecting affine (from rasterio)
  Downloading affine-2.4.0-py3-none-any.whl (15 kB)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2022.12.7)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.3)
Collecting cligj>=0.5 (from rasterio)
  Downloading cligj-0.7.2-py3-none-any.whl (7.1 kB)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.22.4)
Collecting snuggs>=1.4.1 (from rasterio)
  Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)
Collecting click-plugins (from rasterio)
  Downloading click_plugins-1.1.1-py2.py3-none-any.whl (7.5 kB)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (67.7.2)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.0.9)
Installing collected packages: snuggs, cligj, click-plugins, affine, rasterio
Successfully installed affine-2.4.0 click-plugins-1.1.1 cligj-0.7.2 rasterio-1.3.6 snuggs-1.4.7
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting geopandas
  Downloading geopandas-0.13.0-py3-none-any.whl (1.1 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.1/1.1 MB 18.6 MB/s eta 0:00:00
Collecting fiona>=1.8.19 (from geopandas)
  Downloading Fiona-1.9.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (16.5 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 16.5/16.5 MB 36.6 MB/s eta 0:00:00
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from geopandas) (23.1)
Requirement already satisfied: pandas>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.5.3)
Collecting pyproj>=3.0.1 (from geopandas)
  Downloading pyproj-3.5.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (7.7 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.7/7.7 MB 124.7 MB/s eta 0:00:00
Requirement already satisfied: shapely>=1.7.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (2.0.1)
Requirement already satisfied: attrs>=19.2.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2022.12.7)
Requirement already satisfied: click~=8.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (8.1.3)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2022.7.1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (1.22.4)
Installing collected packages: pyproj, fiona, geopandas
Successfully installed fiona-1.9.4 geopandas-0.13.0 pyproj-3.5.0
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting rasterstats
  Downloading rasterstats-0.18.0-py3-none-any.whl (17 kB)
Requirement already satisfied: affine<3.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.4.0)
Requirement already satisfied: click>7.1 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (8.1.3)
Requirement already satisfied: cligj>=0.4 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (0.7.2)
Collecting fiona<1.9 (from rasterstats)
  Downloading Fiona-1.8.22-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (16.6 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 16.6/16.6 MB 117.5 MB/s eta 0:00:00
Requirement already satisfied: numpy>=1.9 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.22.4)
Requirement already satisfied: rasterio>=1.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.3.6)
Collecting simplejson (from rasterstats)
  Downloading simplejson-3.19.1-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (137 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 137.9/137.9 kB 22.3 MB/s eta 0:00:00
Requirement already satisfied: shapely in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.0.1)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (2022.12.7)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.1.1)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.16.0)
Collecting munch (from fiona<1.9->rasterstats)
  Downloading munch-3.0.0-py2.py3-none-any.whl (10 kB)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (67.7.2)
Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from rasterio>=1.0->rasterstats) (1.4.7)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio>=1.0->rasterstats) (3.0.9)
Installing collected packages: simplejson, munch, fiona, rasterstats
  Attempting uninstall: fiona
    Found existing installation: Fiona 1.9.4
    Uninstalling Fiona-1.9.4:
      Successfully uninstalled Fiona-1.9.4
Successfully installed fiona-1.8.22 munch-3.0.0 rasterstats-0.18.0 simplejson-3.19.1
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: earthengine-api in /usr/local/lib/python3.10/dist-packages (0.1.350)
Requirement already satisfied: google-cloud-storage in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.8.0)
Requirement already satisfied: google-api-python-client>=1.12.1 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.84.0)
Requirement already satisfied: google-auth>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.17.3)
Requirement already satisfied: google-auth-httplib2>=0.0.3 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (0.1.0)
Requirement already satisfied: httplib2<1dev,>=0.9.2 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (0.21.0)
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.27.1)
Requirement already satisfied: google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5 in /usr/local/lib/python3.10/dist-packages (from google-api-python-client>=1.12.1->earthengine-api) (2.11.0)
Requirement already satisfied: uritemplate<5,>=3.0.1 in /usr/local/lib/python3.10/dist-packages (from google-api-python-client>=1.12.1->earthengine-api) (4.1.1)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (5.3.0)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (0.3.0)
Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (1.16.0)
Requirement already satisfied: rsa<5,>=3.1.4 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (4.9)
Requirement already satisfied: pyparsing!=3.0.0,!=3.0.1,!=3.0.2,!=3.0.3,<4,>=2.4.2 in /usr/local/lib/python3.10/dist-packages (from httplib2<1dev,>=0.9.2->earthengine-api) (3.0.9)
Requirement already satisfied: google-cloud-core<3.0dev,>=2.3.0 in /usr/local/lib/python3.10/dist-packages (from google-cloud-storage->earthengine-api) (2.3.2)
Requirement already satisfied: google-resumable-media>=2.3.2 in /usr/local/lib/python3.10/dist-packages (from google-cloud-storage->earthengine-api) (2.5.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (3.4)
Requirement already satisfied: googleapis-common-protos<2.0dev,>=1.56.2 in /usr/local/lib/python3.10/dist-packages (from google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5->google-api-python-client>=1.12.1->earthengine-api) (1.59.0)
Requirement already satisfied: protobuf!=3.20.0,!=3.20.1,!=4.21.0,!=4.21.1,!=4.21.2,!=4.21.3,!=4.21.4,!=4.21.5,<5.0.0dev,>=3.19.5 in /usr/local/lib/python3.10/dist-packages (from google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5->google-api-python-client>=1.12.1->earthengine-api) (3.20.3)
Requirement already satisfied: google-crc32c<2.0dev,>=1.0 in /usr/local/lib/python3.10/dist-packages (from google-resumable-media>=2.3.2->google-cloud-storage->earthengine-api) (1.5.0)
Requirement already satisfied: pyasn1<0.6.0,>=0.4.6 in /usr/local/lib/python3.10/dist-packages (from pyasn1-modules>=0.2.1->google-auth>=1.4.1->earthengine-api) (0.5.0)
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (2.27.1)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests) (3.4)
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting sentinelsat
  Downloading sentinelsat-1.2.1-py3-none-any.whl (48 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 48.8/48.8 kB 3.2 MB/s eta 0:00:00
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (2.27.1)
Requirement already satisfied: click>=7.1 in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (8.1.3)
Collecting html2text (from sentinelsat)
  Downloading html2text-2020.1.16-py3-none-any.whl (32 kB)
Collecting geojson>=2 (from sentinelsat)
  Downloading geojson-3.0.1-py3-none-any.whl (15 kB)
Requirement already satisfied: tqdm>=4.58 in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (4.65.0)
Collecting geomet (from sentinelsat)
  Downloading geomet-1.0.0-py3-none-any.whl (28 kB)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from geomet->sentinelsat) (1.16.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (3.4)
Installing collected packages: html2text, geomet, geojson, sentinelsat
Successfully installed geojson-3.0.1 geomet-1.0.0 html2text-2020.1.16 sentinelsat-1.2.1

Setting up some directory paths on Google Drive

A shapefile of Larache province is in my Google Drive. I drew a polygon and saved it as a shapefile on http://www.geojson.io.

# Connecting to Google Earth Engine API
# This opens a web page where i  enter my account information and a verification code is also provided. This code is what i paste into the terminal.
!earthengine authenticate

ee.Initialize()

To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does not start automatically, please manually browse the URL below.

    https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=m6jhmZYUMFpwsVc7HGY_G9JPMbQ8rSlu38E6KuTfon8&tc=qv5LaParODYGIYSxAjbHdPXmk0FuKvtCKGY7nEk4YsE&cc=It22cv2Oj5kPlxLjRdYc79Gu7DVH3ZkLsxjBNxXapOc

The authorization workflow will generate a code, which you should paste in the box below.
Enter verification code: 4/1AbUR2VM1jcMN2HebXQz4IYDnKP9_8Xg8RInpG8gTYo-ZrsMEsR5Pdi9gjSw

Successfully saved authorization token.

'''
--------------------------------------------------------------
The following block of code originates  from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''

# set up your directories for the satellite data
# Note that we do all the downloading and data analysis on the temporary drive
#    on Colab. We will copy the output directory to our Google Drive at the end.
#    Colab has more disk space (about 40 GB free space) than Google Drive (15 GB).
#    However, the data on the Colab disk space are NOT kept when you log out.

# path to your Google Drive
print("Connected to data directory: " + rootdir)

# path to your temporary drive on the Colab Virtual Machine. This will be removed
#   when the session ends.
cd = "/content/work"

# directory for downloading the Sentinel-2 composites
# Note that we are using the 'join' function imported from the os library here
# It is an easy way of merging strings into a directory structure.
# It is clever and chooses the / or \ depending on whether you are on Windows or Linux.
downloaddir = join(cd, 'download') # where we save the downloaded images

# CAREFUL: This code removes the named directories and everything inside them to free up space
# Because the downloaddir is in your temporary drive, it should be empty at the
#    start of the session.
# Note: shutil provides a lot of useful functions for file and directory management
try:
  shutil.rmtree(downloaddir)
except:
  print(downloaddir + " not found.")

# create the new directories, unless they already exist
os.makedirs(cd, exist_ok=True)
os.makedirs(downloaddir, exist_ok=True)

print("Connected to Colab temporary data directory: " + cd)

print("\nList of contents of " + rootdir)
for f in sorted(os.listdir(rootdir)):
  print(f)

Connected to data directory: /content/drive/MyDrive/satelliteCW1
/content/work/download not found.
Connected to Colab temporary data directory: /content/work

List of contents of /content/drive/MyDrive/satelliteCW1
.ipynb_checkpoints
229010645_GY7709_CW1 BACKUP2.ipynb
229010645_GY7709_CW1 BACKUP3.ipynb
229010645_GY7709_CW1.html
229010645_GY7709_CW1.ipynb
229010645_GY7709_CW1_backup.ipynb
229010645_GY7709_CW1_original.html
229010645_GY7709_CW2.html
229010645_GY7709_CW2.ipynb
NBR_larache_after.tiff
NBR_larache_before.tiff
NDVI_larache_after.tiff
NDVI_larache_before.tiff
__pycache__
backup
burned_layers
dNBR.tif
dNBR.tiff
dNBR_warped.tif
larache_after
larache_before
larache_before_fire.tif
larache_before_fire_warped.tif
larache_new_shapefile
larache_shapefiles
leicestershire
oakham
practical_week32_SentinelSat.ipynb
pygge.ipynb
pygge.py
sencredentials.txt
taza
unburned_layers

Defining search parameters for the before fire image of Larache

modifying some of the parameters and uploading the shapefile of Larache called "POLYGON.shp" .

'''
--------------------------------------------------------------
The following block of code  is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# EDITING THE SEARCH OPTIONS BELOW

shapefile = join(rootdir, 'larache_new_shapefile', 'FULL_POLYGON.shp') # ESRI Shapefile of the study area

# checking whether the shapefile exists
if os.path.exists(shapefile):
  print('Shapefile found: '+shapefile)
else:
  print('ERROR: Shapefile not found: '+shapefile)
  print('Upload a shapefile to your Google Drive directory: '+ rootdir)

# Defining a date range for the search
datefrom = '2022-07-01' # start date for imagery search
dateto   = '2022-07-10' # end date for imagery search
time_range = [datefrom, dateto] # format as a list

# Defining which cloud cover percentage we accept in the images
clouds = 8 # maximum acceptable cloud cover in 

Shapefile found: /content/drive/MyDrive/satelliteCW1/larache_new_shapefile/FULL_POLYGON.shp

Getting some information about the shapefile.

'''
--------------------------------------------------------------
The following block of code originates from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Getting the shapefile layer's extent, CRS and EPSG code
extent, outSpatialRef, epsg = pygge.get_shp_extent(shapefile)
print("Extent of the area of interest (shapefile):\n", extent)
print(type(extent))
print("\nCoordinate referencing system (CRS) of the shapefile:\n", outSpatialRef)
print('EPSG code: ', epsg)

Extent of the area of interest (shapefile):
 (-6.08, -6.02, 35.21, 35.31)
<class 'tuple'>

Coordinate referencing system (CRS) of the shapefile:
 GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]
EPSG code:  4326

Getting the extent of the shapefile into a format that Google Earth Engine understands.

'''
--------------------------------------------------------------
The following block of code originates  from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# GEE needs a special format for defining an area of interest. 
# It has to be a GeoJSON Polygon and the coordinates should be first defined in a list and then converted using ee.Geometry. 
extent_list = list(extent)
print(extent_list)
print(type(extent_list))
# close the list of polygon coordinates by adding the starting node at the end again
# and make list elements in the form of coordinate pairs (y,x)
area_list = list([(extent[0], extent[2]),(extent[1], extent[2]),(extent[1], extent[3]),(extent[0], extent[3]),(extent[0], extent[2])])
print(area_list)
print(type(area_list))

search_area = ee.Geometry.Polygon(area_list)
print(search_area)
print(type(search_area))

[-6.08, -6.02, 35.21, 35.31]
<class 'list'>
[(-6.08, 35.21), (-6.02, 35.21), (-6.02, 35.31), (-6.08, 35.31), (-6.08, 35.21)]
<class 'list'>
ee.Geometry({
  "functionInvocationValue": {
    "functionName": "GeometryConstructors.Polygon",
    "arguments": {
      "coordinates": {
        "constantValue": [
          [
            [
              -6.08,
              35.21
            ],
            [
              -6.02,
              35.21
            ],
            [
              -6.02,
              35.31
            ],
            [
              -6.08,
              35.31
            ],
            [
              -6.08,
              35.21
            ]
          ]
        ]
      },
      "evenOdd": {
        "constantValue": true
      }
    }
  }
})
<class 'ee.geometry.Geometry'>

Accessing the Sentinel-2 collection on Google Earth Engine and running our search.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Obtaining download links for image composites from an image collection on Google Earth Engine


# Name of the Sentinel 2 image collection
s2collection = ('COPERNICUS/S2')

# getting the median composite of Sentinel-2 images in the time range
s2median = pygge.obtain_image_sentinel(s2collection, time_range, search_area, clouds)

# Downloading the followig bands as a list of strings namely; Blue, Green,Red,NIR and SWIR

bands = [ 'B2','B3','B4', 'B8','B12']
print(bands)

# spatial resolution of the downloaded data
resolution = 20 # in units of metres

# Download images in Geotiff, using the get_url(name, image, scale, region) method
# ‘region’ is obtained from the area, but the format is adjusted using get_region(geom) method
search_region = pygge.get_region(search_area)
s2url = pygge.get_url('s2', s2median.select(bands), resolution, search_region, filePerBand=False)
print(s2url)

['B2', 'B3', 'B4', 'B8', 'B12']
https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/44e3d97d9adfca34655516a4c2fea849-11e933b4d0f41ebc8d8b877d0589bca7:getPixels

Download the data for the before fire image of Larache

The next cell downloads the image composite as a zip file and unzips it.

'''
--------------------------------------------------------------
The following block of code originates  from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# change directory to download directory
os.chdir(downloaddir)

# requesting information on the file to be downloaded
f = pygge.requests.get(s2url, stream =True)

# checking whether it is a zip file
check = zipfile.is_zipfile(io.BytesIO(f.content))

# either download the file as is, or unzip it
while not check:
    f = requests.get(s2url, stream =True)
    check = zipfile.is_zipfile(io.BytesIO(f.content))
else:
    z = zipfile.ZipFile(io.BytesIO(f.content))
    z.extractall()

Exploring the data directory structure of the downloaded files

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# where  the downloaded Sentinel-2 images are stored
os.chdir(downloaddir)
print("contents of ", downloaddir, ":")
!ls -l

contents of  /content/work/download :
total 2120
-rw-r--r-- 1 root root 2168674 May 12 13:28 s2.tif

The downloaded file " s2.tif " is seen.

the downloaded images are saved to a temporary directory that will be deleted when the virtual machine is closed. To save the images to my local directory, this is how it went;

I went to my Google Colab folder in the panel on the left hand side.

Found the download directory and clicked on a Sentinel-2 image folder.

Right-clicked on it ,renamed it and selected 'download' to save it.

Showing the before fire image of Larache as a colour composite

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# getting list of all tiff files in the directory
allfiles = [f for f in listdir(downloaddir) if isfile(join(downloaddir, f))]
print(allfiles)

# selecting the file for visualisation
thisfile = allfiles[0]
print(thisfile)

['s2.tif']
s2.tif

True Colour Composite of before fire image of Larache

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure 
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')

# the downloaded file is float32 data format
# for plotting,  uint8 data format is needed

# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[3,2,1], percentiles=[0,99])

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

False Colour Composite of before fire image of Larache

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')

# the downloaded file is float32 data format
# for plotting,  uint8 data format is needed

# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[4,2,1], percentiles=[0,99])

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

Warping the downloaded before fire image composite of Larache into another map projection

The coordinate reference system (CRS) of the downloaded image composite is not in the UK national map projection. hence it is reprojected.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing the EPSG code of our shapefile into which we want to reproject the TCI images
print("Reprojecting image to EPSG projection ", epsg)

# making a file name for the new file
warpfile = thisfile.split(sep='.')[0] + '_warped.tif'
print("We are in this directory: ")
!pwd
print("Input file: ", thisfile)
print("Output file: ", warpfile)

# calling the easy_warp function
tmp = pygge.easy_warp(thisfile, warpfile, epsg)

Reprojecting image to EPSG projection  4326
We are in this directory: 
/content/work/download

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

Input file:  s2.tif
Output file:  s2_warped.tif
Creating warped file:s2_warped.tif

Plotting the shapefile on top of the raster

To see the locations of our polygons on top of our image composite, we can do that with the Geopandas library.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with subplots
fig, ax = plt.subplots(2,1, figsize=(10,16))
fig.patch.set_facecolor('white')

# plotting the image with full extent in true colour
pygge.easy_plot(warpfile, ax=ax[0], percentiles=[0,98], bands=[3,2,1],
                shapefile=shapefile, fillcolor="none", linecolor="yellow", 
                title="LARACHE BEFORE FIRE")

# plotting the image with full extent in false colour
pygge.easy_plot(warpfile, ax=ax[1], percentiles=[0,98], bands=[4,2,1],
                shapefile=shapefile, fillcolor="none", linecolor="yellow", 
                title="LARACHE BEFORE FIRE")

Defining the search parameters for the after fire image of Larache

We change the search dates to get a new image after the fire had burnt out. Fires were first reported on 15 July 2022

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# EDITing THE SEARCH OPTIONS BELOW

shapefile = join(rootdir, 'larache_shapefiles', 'FULL_POLYGON.shp') # ESRI Shapefile of the study area

# checking whether the shapefile exists
if os.path.exists(shapefile):
  print('Shapefile found: '+shapefile)
else:
  print('ERROR: Shapefile not found: '+shapefile)
  print('Upload a shapefile to your Google Drive directory: '+ rootdir)

# Defining a date range for the search
datefrom = '2022-07-30' # start date for imagery search
dateto   = '2022-08-05' # end date for imagery search
time_range = [datefrom, dateto] # format as a list

# Defining percentage cloud cover accepted in the images
clouds = 10 # maximum acceptable cloud cover in %

Shapefile found: /content/drive/MyDrive/satelliteCW1/larache_shapefiles/FULL_POLYGON.shp

Get some information about our shapefile.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Getting the shapefile layer's extent, CRS and EPSG code
extent, outSpatialRef, epsg = pygge.get_shp_extent(shapefile)
print("Extent of the area of interest (shapefile):\n", extent)
print(type(extent))
print("\nCoordinate referencing system (CRS) of the shapefile:\n", outSpatialRef)
print('EPSG code: ', epsg)

Extent of the area of interest (shapefile):
 (-6.08, -6.02, 35.21, 35.31)
<class 'tuple'>

Coordinate referencing system (CRS) of the shapefile:
 GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]
EPSG code:  4326

Getting the extent of the shapefile into a format that Google Earth Engine understands.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# GEE needs a special format for defining an area of interest. 
# It has to be a GeoJSON Polygon and the coordinates should be first defined in a list and then converted using ee.Geometry. 
extent_list = list(extent)
print(extent_list)
print(type(extent_list))
# close the list of polygon coordinates by adding the starting node at the end again
# and make list elements in the form of coordinate pairs (y,x)
area_list = list([(extent[0], extent[2]),(extent[1], extent[2]),(extent[1], extent[3]),(extent[0], extent[3]),(extent[0], extent[2])])
print(area_list)
print(type(area_list))

search_area = ee.Geometry.Polygon(area_list)
print(search_area)
print(type(search_area))

[-6.08, -6.02, 35.21, 35.31]
<class 'list'>
[(-6.08, 35.21), (-6.02, 35.21), (-6.02, 35.31), (-6.08, 35.31), (-6.08, 35.21)]
<class 'list'>
ee.Geometry({
  "functionInvocationValue": {
    "functionName": "GeometryConstructors.Polygon",
    "arguments": {
      "coordinates": {
        "constantValue": [
          [
            [
              -6.08,
              35.21
            ],
            [
              -6.02,
              35.21
            ],
            [
              -6.02,
              35.31
            ],
            [
              -6.08,
              35.31
            ],
            [
              -6.08,
              35.21
            ]
          ]
        ]
      },
      "evenOdd": {
        "constantValue": true
      }
    }
  }
})
<class 'ee.geometry.Geometry'>

Accessing the Sentinel-2 collection on Google Earth Engine and running the search.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Obtain download links for image composites from an image collection on Google Earth Engine

# Name of the Sentinel 2 image collection
s2collection = ('COPERNICUS/S2')

# getting the median composite of Sentinel-2 images in the time range
s2median = pygge.obtain_image_sentinel(s2collection, time_range, search_area, clouds)

# Blue, Green,Red,NIR and SWIR are downloaded
bands = ['B2','B3','B4', 'B8','B12' ]
print(bands)

# spatial resolution of the downloaded data
resolution = 20 # in units of metres

# Downloading images in Geotiff, using the get_url(name, image, scale, region) method
# ‘region’ is obtained from the area, but the format has to be adjusted using get_region(geom) method
search_region = pygge.get_region(search_area)
s2url = pygge.get_url('s2', s2median.select(bands), resolution, search_region, filePerBand=False)
print(s2url)

['B2', 'B3', 'B4', 'B8', 'B12']
https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/3b62a644518dacef2cf989403d8c3e3f-44432a28d67b784f448895200db8318e:getPixels

Downloading the data

The next cell downloads the image composite as a zip file and unzips it.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# changing directory to download directory
os.chdir(downloaddir)

# requesting information on the file to be downloaded
f = pygge.requests.get(s2url, stream =True)

# checking whether it is a zip file
check = zipfile.is_zipfile(io.BytesIO(f.content))

# either download the file as is, or unzip it
while not check:
    f = requests.get(s2url, stream =True)
    check = zipfile.is_zipfile(io.BytesIO(f.content))
else:
    z = zipfile.ZipFile(io.BytesIO(f.content))
    z.extractall()

Exploring the data directory structure of our downloaded files

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# where we stored the downloaded Sentinel-2 images
os.chdir(downloaddir)
print("contents of ", downloaddir, ":")
!ls -l

contents of  /content/work/download :
total 2200
-rw-r--r-- 1 root root 2251399 Apr 27 11:10 s2.tif

the downloaded file is " s2.tif ".

the downloaded images are saved to a temporary directory that will be deleted when the virtual machine is closed. To save the images to my local directory, this is how it went;

I went to my Google Colab folder in the panel on the left hand side.

Found the download directory and clicked on a Sentinel-2 image folder.

Right-clicked on it ,renamed it and selected 'download' to save it.

Showing the after fire image of Larache as a colour composite

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# getting list of all tiff files in the directory
allfiles = [f for f in listdir(downloaddir) if isfile(join(downloaddir, f))]
print(allfiles)

# selecting the file for visualisation
thisfile = allfiles[0]
print(thisfile)

['s2.tif']
s2.tif

True Colour composite of after fire image of Larache

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')

# the downloaded file is float32 data format
# for plotting,  uint8 data format is needed

# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[3,2,1], percentiles=[0,99])

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

False Colour Composite of after fire image of Larache

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')

# the downloaded file is float32 data format
# for plotting, we need uint8 data format

# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[4,2,1], percentiles=[0,99])

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

Warping the downloaded after fire image composite of larache into another map projection

The coordinate reference system (CRS) of the downloaded image composite is not in the UK national map projection. It is reprojected.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing the EPSG code of our shapefile into which we want to reproject the TCI images
print("Reprojecting image to EPSG projection ", epsg)

# making a file name for our new file
warpfile = thisfile.split(sep='.')[0] + '_warped.tif'
print("We are in this directory: ")
!pwd
print("Input file: ", thisfile)
print("Output file: ", warpfile)

# calling the easy_warp function
tmp = pygge.easy_warp(thisfile, warpfile, epsg)

Reprojecting image to EPSG projection  4326
We are in this directory: 
/content/work/download

WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

Input file:  s2.tif
Output file:  s2_warped.tif
Creating warped file:s2_warped.tif

Plotting the shapefile on top of the raster

to see the locations of our polygons on top of our image composite,We do that with the Geopandas library.

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with subplots
fig, ax = plt.subplots(2,1, figsize=(10,16))
fig.patch.set_facecolor('white')

# plotting the image with full extent
pygge.easy_plot(warpfile, ax=ax[0], percentiles=[0,98], bands=[3,2,1],
                shapefile=shapefile, fillcolor="none", linecolor="yellow", 
                title="LARACHE AFTER FIRE")

# plotting the image with full extent
pygge.easy_plot(warpfile, ax=ax[1], percentiles=[0,98], bands=[4,2,1],
                shapefile=shapefile, fillcolor="none", linecolor="yellow", 
                title="LARACHE AFTER FIRE")



                               

READING IN DOWNLOADED TIF FILE OF LARACHE BEFORE FIRE AND GETTING ITS METADATA

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# reading in the warped before fire downloaded image of larache and checking its cordinate reference system
larache_before = rasterio.open(join(rootdir +  "/larache_before/before_fire_warped.tif"))
src1=larache_before
d1=src1.crs.to_wkt()
print(type(d1))
pprint(d1)

<class 'str'>
('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS '
 '84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')

READING IN DOWNLOADED TIF FILE OF LARACHE AFTER FIRE AND GETTING ITS METADATA

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# reading in the warped after fire downloaded image of larache and checking its cordinate system
larache_after = rasterio.open(join(rootdir +  "/larache_after/after_fire_warped.tif"))
src2=larache_after
d2=src2.crs.to_wkt()
print(type(d2))
pprint(d2)

<class 'str'>
('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS '
 '84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')

CHECKING IF THE BEFORE AND AFTER DATASETS ARE OF SAME SHAPE

# checking if the before and after image data are the same shape?
larache_before.shape == larache_after.shape

True

CALCULATION OF NBR FOR THE BEFORE FIRE IMAGE OF LARACHE

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
swir_before = larache_before.read(5) # band in position 5 in our stacked image i.e SWIR
nir_before = larache_before.read(4) # band in position 4 in our stacked image i.e NIR

# calculation of NBR before the fire as flaoting point array
nbr_before = (nir_before.astype(float) - swir_before.astype(float)) / (nir_before.astype(float) + swir_before.astype(float))


# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# printing some image statistics, ignoring missing values (nan)
print("minimum NBR_before = ", np.nanmin(nbr_before))
print("mean NBR_before = ", np.nanmean(nbr_before))
print("maximum NBR_before = ", np.nanmax(nbr_before))
print("standard deviation = ", np.nanstd(nbr_before))

minimum NBR_before =  -0.14298390392538413
mean NBR_before =  0.10239791290152904
maximum NBR_before =  0.44895833342572644
standard deviation =  0.0735262823032699

PLOT OF NBR BEFORE THE FIRE

plt.imshow(nbr_before, cmap='PiYG')    # colours to use  include Purple and Green
plt.colorbar()
plt.title('NBR of Larache before fire')  # title of plot
plt.show()                               # plot

CALCULATION OF NBR FOR THE AFTER FIRE IMAGE OF LARACHE

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
swir_after = larache_after.read(5) # band in position 5 in our stacked image i.e SWIR
nir_after = larache_after.read(4) # band in position 4 in our stacked image i.e NIR

# calculation of NBR after the fire as flaoting point array
nbr_after = (nir_after.astype(float) - swir_after.astype(float)) / (nir_after.astype(float) + swir_after.astype(float))

# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# printing some image statistics, ignoring missing values (nan)
print("minimum NBR_after = ", np.nanmin(nbr_after))
print("mean NBR_after = ", np.nanmean(nbr_after))
print("maximum NBR_after = ", np.nanmax(nbr_after))
print("standard deviation = ", np.nanstd(nbr_after))

minimum NBR_after =  -0.3224225354386313
mean NBR_after =  0.049028595859026615
maximum NBR_after =  0.5460599474035327
standard deviation =  0.11624440093712922

PLOT OF NBR AFTER THE FIRE

plt.imshow(nbr_after, cmap='PiYG') # Purple and Green selected as colours in plot
plt.colorbar()
plt.title('NBR of Larache after fire')  # title of plot
plt.show()                              # display plot or map

CALCULATION OF dNBR AS DIFFERENCE OF BEFORE AND AFTER NBR RASTERS

# finding the difference between the NBR before and after the fire
dNBR = nbr_before - nbr_after

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''

# printing some image statistics, ignoring missing values (nan)
print("minimum dNBR = ", np.nanmin(dNBR))
print("mean dNBR = ", np.nanmean(dNBR))
print("maximum dNBR = ", np.nanmax(dNBR))
print("standard deviation = ", np.nanstd(dNBR))

minimum dNBR =  -0.363571844669598
mean dNBR =  0.05336931704250243
maximum dNBR =  0.5638484811508442
standard deviation =  0.10752007563079775

MAP SHOWING AREAS WHERE THE FIRE DAMAGE HAS BEEN GREATEST

# Defining the dNBR values and the corresponding severity levels
dNBR

levels = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
labels = ['Enhanced regrowth,high','Enhanced regrowth,low','Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity']

# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['purple','cyan','green', 'yellow', 'blue', 'darkred', 'red'])

# Setting the boundaries and norm of the color map
bounds = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
norm = colors.BoundaryNorm(bounds, cmap.N)

# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))

# Plotting the dNBR using the custom color map and norm
im = ax.imshow(dNBR, cmap=cmap, norm=norm)

# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)

# Setting the title and axis labels
ax.set_title('dNBR Map with Severity Levels')

plt.show()

CALCULATION OF AREAS OF BURNED CLASSES USING dNBR

# Define the area of each pixel in hectares assuming a spatial resolution of 20 meters
pixel_area = 0.0004 

# Calculate the number of pixels in each class or label
regrowth_high_pixels = np.sum((dNBR >= -0.5) & (dNBR < -0.251))
regrowth_low_pixels = np.sum((dNBR >= -0.25) & (dNBR < -0.101))
unburned_pixels = np.sum((dNBR >= -0.100) & (dNBR < 0.99))
low_severity_pixels = np.sum((dNBR >= 0.100) & (dNBR < 0.269))
mod_low_severity_pixels = np.sum((dNBR >= 0.270) & (dNBR < 0.439))
mod_high_severity_pixels = np.sum((dNBR >= 0.44) & (dNBR < 0.659))
high_severity_pixels = np.sum((dNBR >= 0.660) & (dNBR < 1.300))

# Calculate the area of each class or label in hectares
regrowth_high_area = regrowth_high_pixels * pixel_area
regrowth_low_area = regrowth_low_pixels * pixel_area
unburned_area = unburned_pixels * pixel_area
low_severity_area = low_severity_pixels * pixel_area
mod_low_severity_area = mod_low_severity_pixels * pixel_area
mod_high_severity_area = mod_high_severity_pixels * pixel_area
high_severity_area = high_severity_pixels * pixel_area

# Print the area of each class or label
print('Enhanced regrowth,high area:', regrowth_high_area, 'hectares')
print('Enhanced regrowth, low area:', regrowth_low_area, 'hectares')
print('Unburned area:', unburned_area, 'hectares')
print('Low severity area:', low_severity_area, 'hectares')
print('Moderate-low severity area:', mod_low_severity_area, 'hectares')
print('Moderate-high severity area:', mod_high_severity_area, 'hectares')
print('High severity area:', high_severity_area, 'hectares')

Enhanced regrowth,high area: 0.0036000000000000003 hectares
Enhanced regrowth, low area: 0.0664 hectares
Unburned area: 74.6992 hectares
Low severity area: 2.0148 hectares
Moderate-low severity area: 4.2816 hectares
Moderate-high severity area: 1.4420000000000002 hectares
High severity area: 0.0 hectares

CALCULATION OF AVERAGE NBR2 FOR BURNED POLYGON

# Opening the TIF file and the shapefile of burned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_after/after_fire_warped.tif') as src:
    shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/burned_layers/POLYGON.shp')

    # Masking the TIF file using the shapefile
    masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)

    # Calculating the NBR2_burned
    # the NIR band has index 3 and SWIR band has index 4
    NBR2_burned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])

    # Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# print average NBR2 for burned poloygon, ignoring missing values (nan)

print("average NBR2_burned = ", np.nanmean(NBR2_burned))

average NBR2_burned =  -0.13870606

CALCULATION OF AVERAGE NBR1 FOR BURNED POLYGON

# Opening the TIF file and the shapefile of burned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_before/before_fire_warped.tif') as src:
    shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/burned_layers/POLYGON.shp')

    # Masking the TIF file using the shapefile
    masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)

    # Calculating the NBR1_burned
    # the NIR  band has index 3 and SWIR band has index 4
    NBR1_burned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])

    # Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# printing average NBR1 for burned polygon, ignoring missing values (nan)

print("average NBR1_burned = ", np.nanmean(NBR1_burned))

average NBR1_burned =  0.1715795

CALCULATION OF dNBR FOR BURNED POLYGON

#  finding the difference between the NBR before and after the fire
dNBR_burned = NBR1_burned - NBR2_burned

# printing average dNBR for burned polygon, ignoring missing values (nan)

print("average dNBR_burned = ", np.nanmean(dNBR_burned))

average dNBR_burned =  0.3102856

CALCULATION OF AVERAGE NBR1 FOR UNBURNED POLYGON

# Opening the TIF file and the shapefile of unburned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_before/before_fire_warped.tif') as src:
    shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/unburned_layers/POLYGON.shp')

    # Masking the TIF file using the shapefile
    masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)

    # Calculating the NBR1_unburned
    # the NIR band has index 3 and SWIR band has index 4
    NBR1_unburned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])

    # Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# printing average NBR1 for unburned poloygon, ignoring missing values (nan)

print("average NBR1_unburned = ", np.nanmean(NBR1_unburned))

average NBR1_unburned =  0.17309634

CALCULATION OF AVERAGE NBR2 FOR UNBURNED POLYGON

# Opening the TIF file and the shapefile of unburned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_after/after_fire_warped.tif') as src:
    shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/unburned_layers/POLYGON.shp')

    # Masking the TIF file using the shapefile
    masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)

    # Calculating the NBR2_unburned
    # the NIR band has index 3 and SWIR band has index 4
    NBR2_unburned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])

    # Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# printing average NBR2 for unburned poloygon, ignoring missing values (nan)

print("average NBR2_unburned = ", np.nanmean(NBR2_unburned))

average NBR2_unburned =  0.17163754

CALCULATION OF dNBR FOR UNBURNED POLYGON

# finding the difference between the NBR before and after the fire
dNBR_unburned = NBR1_unburned - NBR2_unburned

# printing average dNBR for burned polygon, ignoring missing values (nan)

print("average dNBR_unburned = ", np.nanmean(dNBR_unburned))

average dNBR_unburned =  0.0014588113

CALCULATION OF NDVI FOR THE BEFORE FIRE IMAGE OF LARACHE

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
r_before = larache_before.read(3) # band 3 in our stacked image is Red
nir_before = larache_before.read(4) # band 4 in our stacked image is NIR

# calculation of NDVI before the fire as flaoting point array
NDVI_before = (nir_before.astype(float) - r_before.astype(float)) / (nir_before.astype(float) + r_before.astype(float))

# Ignore division by zero
np.seterr(divide='ignore', invalid='ignore')

# print some image statistics, ignoring missing values (nan)
print("minimum NDVI_before = ", np.nanmin(NDVI_before))
print("mean NDVI_before = ", np.nanmean(NDVI_before))
print("maximum NDVI_before = ", np.nanmax(NDVI_before))
print("standard deviation = ", np.nanstd(NDVI_before))

minimum NDVI_before =  -0.0031236251459653844
mean NDVI_before =  0.16960801456931313
maximum NDVI_before =  0.44798869776436046
standard deviation =  0.053921411038247144

CALCULATION OF NDVI FOR THE AFTER FIRE IMAGE OF LARACHE

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
r_after = larache_after.read(3) # band 3 in our stacked image is Red
nir_after = larache_after.read(4) # band 4 in our stacked image is NIR

# calculation of NDVI after the fire as flaoting point array
NDVI_after = (nir_after.astype(float) - r_after.astype(float)) / (nir_after.astype(float) + r_after.astype(float))

# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')

# print some image statistics, ignoring missing values (nan)
print("minimum NDVI_after = ", np.nanmin(NDVI_after))
print("mean NDVI_after = ", np.nanmean(NDVI_after))
print("maximum NDVI_after = ", np.nanmax(NDVI_after))
print("standard deviation = ", np.nanstd(NDVI_after))

minimum NDVI_after =  -0.06513785137467486
mean NDVI_after =  0.17277971166410844
maximum NDVI_after =  0.5750123655099424
standard deviation =  0.07998356629775627

MAPS OF NDVI BEFORE THE FIRE

plt.imshow(NDVI_before, cmap='RdYlGn') # selection of Red to Green colours
plt.colorbar()                           # colour bar
plt.title('NDVI of Larache before fire')  # map title
plt.show()                                # showing plot

# Defining the pre-fire NDVI  levels
NDVI_before

levels = [-1, 0, 0.2, 0.5]
labels = ['barren or non-vegetated','sparse or stressed vegetation','moderate or healthy vegetation', 'dense or very healthy vegetation']

# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['orange', 'yellow', 'lightgreen', 'darkgreen'])

# Setting the boundaries and norm of the color map
bounds = [-1, 0, 0.2, 0.5]
norm = colors.BoundaryNorm(bounds, cmap.N)

# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))

# Plotting the pre-fire NDVI using the custom color map and norm
im = ax.imshow(NDVI_before, cmap=cmap, norm=norm)

# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)

# Setting the title and axis labels
ax.set_title('Pre-fire NDVI Map of Larache ')

plt.show()

MAPS OF NDVI AFTER THE FIRE

plt.imshow(NDVI_after, cmap='RdYlGn')      # Red to Green colour selected
plt.colorbar()
plt.title('NDVI of Larache after fire')    # title
plt.show()                                  

# Defining the post-fire NDVI  levels
NDVI_after

levels = [-1, 0, 0.2, 0.5]
labels = ['barren or non-vegetated','sparse or stressed vegetation','moderate or healthy vegetation', 'dense or very healthy vegetation']

# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['orange', 'yellow', 'lightgreen', 'darkgreen'])

# Setting the boundaries and norm of the color map
bounds = [-1, 0, 0.2, 0.5]
norm = colors.BoundaryNorm(bounds, cmap.N)

# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))

# Plotting the post-fire NDVI using the custom color map and norm
im = ax.imshow(NDVI_after, cmap=cmap, norm=norm)

# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)

# Setting the title and axis labels
ax.set_title('Post-fire NDVI Map of Larache ')

plt.show()

CALCULATING THE NDVI DIFFERENCE (dNDVI)

# finding the difference between the NDVI before and after the fire
dNDVI = NDVI_before - NDVI_after

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''

# printing some image statistics, ignoring missing values (nan)
print("minimum dNDVI = ", np.nanmin(dNDVI))
print("mean dNDVI = ", np.nanmean(dNDVI))
print("maximum dNDVI = ", np.nanmax(dNDVI))
print("standard deviation = ", np.nanstd(dNDVI))

minimum dNDVI =  -0.3098725340997325
mean dNDVI =  -0.0031716970947953077
maximum dNDVI =  0.33930637628747506
standard deviation =  0.06478898685069448

dNDVI based BURN SEVERITY MAP

# Defining the dNDVI values and the corresponding severity levels
dNDVI

levels = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
labels = ['Enhanced regrowth,high','Enhanced regrowth,low','Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity']

# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['purple','cyan','green', 'yellow', 'blue', 'darkred', 'red'])

# Setting the boundaries and norm of the color map
bounds = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
norm = colors.BoundaryNorm(bounds, cmap.N)

# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))

# Plotting the dNDVI using the custom color map and norm
im = ax.imshow(dNDVI, cmap=cmap, norm=norm)

# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)

# Setting the title and axis labels
ax.set_title('dNDVI Map with Severity Levels')

plt.show()

CALCULATION OF BURNED AREAS BASED ON dNDVI

# Define the area of each pixel in hectares assuming a spatial resolution of 20 meters
pixel_area = 0.0004 

# Calculate the number of pixels in each class or label
regrowth_high_pixels = np.sum((dNDVI >= -0.5) & (dNDVI < -0.251))
regrowth_low_pixels = np.sum((dNDVI >= -0.25) & (dNDVI < -0.101))
unburned_pixels = np.sum((dNDVI >= -0.100) & (dNDVI < 0.99))
low_severity_pixels = np.sum((dNDVI >= 0.100) & (dNDVI < 0.269))
mod_low_severity_pixels = np.sum((dNDVI >= 0.270) & (dNDVI < 0.439))
mod_high_severity_pixels = np.sum((dNDVI >= 0.44) & (dNDVI < 0.659))
high_severity_pixels = np.sum((dNDVI >= 0.660) & (dNDVI < 1.300))

# Calculate the area of each class or label in hectares
regrowth_high_area = regrowth_high_pixels * pixel_area
regrowth_low_area = regrowth_low_pixels * pixel_area
unburned_area = unburned_pixels * pixel_area
low_severity_area = low_severity_pixels * pixel_area
mod_low_severity_area = mod_low_severity_pixels * pixel_area
mod_high_severity_area = mod_high_severity_pixels * pixel_area
high_severity_area = high_severity_pixels * pixel_area

# Print the area of each class or label
print('Enhanced regrowth,high area:', regrowth_high_area, 'hectares')
print('Enhanced regrowth, low area:', regrowth_low_area, 'hectares')
print('Unburned area:', unburned_area, 'hectares')
print('Low severity area:', low_severity_area, 'hectares')
print('Moderate-low severity area:', mod_low_severity_area, 'hectares')
print('Moderate-high severity area:', mod_high_severity_area, 'hectares')
print('High severity area:', high_severity_area, 'hectares')

Enhanced regrowth,high area: 0.004 hectares
Enhanced regrowth, low area: 0.2932 hectares
Unburned area: 74.45960000000001 hectares
Low severity area: 6.163600000000001 hectares
Moderate-low severity area: 0.15760000000000002 hectares
Moderate-high severity area: 0.0 hectares
High severity area: 0.0 hectares

TESTING WHETHER THE dNBR WAS CORRELATED TO dNDVI

import numpy as np
from scipy.stats import linregress
import matplotlib.pyplot as plt

# Assuming dNBR and dNDVI are 2D arrays with the same shape
dNBR_flat = -dNBR.flatten()
dNDVI_flat = dNDVI.flatten()

# Calculate the linear regression line
slope, intercept, r_value, p_value, std_err = linregress(dNBR_flat, dNDVI_flat)
line = slope * dNBR_flat + intercept

# Create a scatter plot
fig, ax = plt.subplots()
ax.scatter(dNBR_flat, dNDVI_flat, alpha=0.5)

# Add the linear regression line to the scatter plot with the equation and correlation coefficient
eqn_label = f'dNDVI = {slope:.2f}(dNBR) + {intercept:.2f}'
corr_label = f'r = {r_value:.2f}'
ax.plot(dNBR_flat, line, color='red', label=f'{eqn_label}\n{corr_label}')

# Set the x-axis and y-axis labels
ax.set_xlabel('dNBR')
ax.set_ylabel('dNDVI')

# Set the title of the plot
ax.set_title('Plot of dNDVI versus dNBR')

# Add the legend to the plot
ax.legend()

plt.show()

'''
--------------------------------------------------------------
The following block of code originates and is modified from:

Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
#To save a notebook as an html file, i add a new code cell at the end of the notebook with the following contents:
!pip install -U notebook-as-pdf
!pyppeteer-install
!jupyter nbconvert /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.ipynb --to html
 

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting notebook-as-pdf
  Downloading notebook_as_pdf-0.5.0-py3-none-any.whl (6.5 kB)
Requirement already satisfied: nbconvert in /usr/local/lib/python3.10/dist-packages (from notebook-as-pdf) (6.5.4)
Collecting pyppeteer (from notebook-as-pdf)
  Downloading pyppeteer-1.0.2-py3-none-any.whl (83 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 83.4/83.4 kB 3.7 MB/s eta 0:00:00
Collecting PyPDF2 (from notebook-as-pdf)
  Downloading pypdf2-3.0.1-py3-none-any.whl (232 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 232.6/232.6 kB 11.9 MB/s eta 0:00:00
Requirement already satisfied: lxml in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (4.9.2)
Requirement already satisfied: beautifulsoup4 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (4.11.2)
Requirement already satisfied: bleach in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (6.0.0)
Requirement already satisfied: defusedxml in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.7.1)
Requirement already satisfied: entrypoints>=0.2.2 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.4)
Requirement already satisfied: jinja2>=3.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (3.1.2)
Requirement already satisfied: jupyter-core>=4.7 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.3.0)
Requirement already satisfied: jupyterlab-pygments in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.2.2)
Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (2.1.2)
Requirement already satisfied: mistune<2,>=0.8.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.8.4)
Requirement already satisfied: nbclient>=0.5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.7.4)
Requirement already satisfied: nbformat>=5.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.8.0)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (23.1)
Requirement already satisfied: pandocfilters>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (1.5.0)
Requirement already satisfied: pygments>=2.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (2.14.0)
Requirement already satisfied: tinycss2 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (1.2.1)
Requirement already satisfied: traitlets>=5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.7.1)
Requirement already satisfied: appdirs<2.0.0,>=1.4.3 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (1.4.4)
Requirement already satisfied: certifi>=2021 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (2022.12.7)
Collecting importlib-metadata>=1.4 (from pyppeteer->notebook-as-pdf)
  Downloading importlib_metadata-6.6.0-py3-none-any.whl (22 kB)
Collecting pyee<9.0.0,>=8.1.0 (from pyppeteer->notebook-as-pdf)
  Downloading pyee-8.2.2-py2.py3-none-any.whl (12 kB)
Requirement already satisfied: tqdm<5.0.0,>=4.42.1 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (4.65.0)
Requirement already satisfied: urllib3<2.0.0,>=1.25.8 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (1.26.15)
Collecting websockets<11.0,>=10.0 (from pyppeteer->notebook-as-pdf)
  Downloading websockets-10.4-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (106 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 106.8/106.8 kB 14.4 MB/s eta 0:00:00
Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.10/dist-packages (from importlib-metadata>=1.4->pyppeteer->notebook-as-pdf) (3.15.0)
Requirement already satisfied: platformdirs>=2.5 in /usr/local/lib/python3.10/dist-packages (from jupyter-core>=4.7->nbconvert->notebook-as-pdf) (3.3.0)
Requirement already satisfied: jupyter-client>=6.1.12 in /usr/local/lib/python3.10/dist-packages (from nbclient>=0.5.0->nbconvert->notebook-as-pdf) (6.1.12)
Requirement already satisfied: fastjsonschema in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert->notebook-as-pdf) (2.16.3)
Requirement already satisfied: jsonschema>=2.6 in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert->notebook-as-pdf) (4.3.3)
Requirement already satisfied: soupsieve>1.2 in /usr/local/lib/python3.10/dist-packages (from beautifulsoup4->nbconvert->notebook-as-pdf) (2.4.1)
Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert->notebook-as-pdf) (1.16.0)
Requirement already satisfied: webencodings in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert->notebook-as-pdf) (0.5.1)
Requirement already satisfied: attrs>=17.4.0 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert->notebook-as-pdf) (23.1.0)
Requirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert->notebook-as-pdf) (0.19.3)
Requirement already satisfied: pyzmq>=13 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (23.2.1)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (2.8.2)
Requirement already satisfied: tornado>=4.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (6.3.1)
Installing collected packages: pyee, websockets, PyPDF2, importlib-metadata, pyppeteer, notebook-as-pdf
Successfully installed PyPDF2-3.0.1 importlib-metadata-6.6.0 notebook-as-pdf-0.5.0 pyee-8.2.2 pyppeteer-1.0.2 websockets-10.4
[INFO] Starting Chromium download.
100% 109M/109M [00:00<00:00, 234Mb/s] 
[INFO] Beginning extraction
[INFO] Chromium extracted to: /root/.local/share/pyppeteer/local-chromium/588429
[NbConvertApp] Converting notebook /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.ipynb to html
[NbConvertApp] Writing 6741948 bytes to /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.html